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(54) Abstract Title: Microseismic determination of location and origin time of a fracture generated by fracturing 
operation in a hydrocarbon well 

(57) In a method of monitoring a hydraulic fracturing 
operation in a hydrocarbon well 11, microseismic 
signals 14 originating from an induced fracture 14 
are recorded by three component geophone arrays 
121,131 in neighbouring wells 12,13. A wavefield 
inversion algorithm is used to determine the 
mechanism, origin time and location of the seismic 
source 14. The algorithm does not require the 
signals to be resolved into P-wave and S-wave data 
before inversion. The algorithm evaluates Green's 
functions, uses an existing velocity model and 
decomposes a moment tensor to yield parameters 
characteristic of the fracture. The recorded signals 
may be bandpass limited to 0-50 Hz or 0-100 Hz. 
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This print takes account of replacement documents submitted after the date of filing to enable the application to comply 
with the formal requirements of the Patents Rules 1995 

Original Printed on Recycled Paper 
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Method for Monitoring Hydraulic Stimulation 

This invention relates to methods for acquiring seismic data 
monitoring hydraulic stimulation such as fracturing rock 
5 . layers to improve hydrocarbon production of a well . More 
specifically it relates to such methods using seismic 
methods to detect the location of an induced fracture. 

BACKGROUND OF THE INVENTION 
10 In certain situations, workers in the oil and gas industry 
perform a procedure known as "hydraulic fracturing" during 
oil exploration and drilling operations. For example, in 
formations where oil or gas cannot be easily or economically 
extracted from the earth, a hydraulic fracturing operation 
15 is commonly performed. Such a hydraulic fracturing operation 
includes pumping in large amounts of fluid to induce cracks 
in the earth, thereby creating pathways via which the oil 
and gas may flow. After a crack is generated, sand or some 
other material is commonly added to the crack, so that when 
20 the earth closes back up after the pressure is released, the 
sand helps to keep the earth parted. The sand then provides 
a conductive pathway for the oil and gas to flow from the 
newly formed fracture 

25 However, the hydraulic fracturing process does not always 

work very well. The reasons for this are relatively unknown. 
In addition, the hydraulic fractures cannot be readily 
observed, since they are typically thousands of feet below 
the surface of the earth. Therefore, members of the oil and 

30 gas industry have sought diagnostic methods to tell where 
the fractures are, how big the fractures are, how far they 
go and how high they grow. Thus, a diagnostic apparatus and 
method for measuring the hydraulic fracture and the rock 
deformation around the fracture are needed. 
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In previous attempts to solve this problem, certain methods 
have been developed for mapping fractures. For example, one 
of these methods involves seismic sensing. In such a seismic 
sensing operation, micro-earthquakes generated by the 
5 fracturing are analyzed by seismic meters, for example, 
accelerometers . 

A recent study on the use of microseismic imaging for 
. fracture stimulation was published by J. T. Rutledge and W. 

10 S. Phillips. In an typical operational setting as described 
in greater detail in FIG. 1 below, three- component geophones 
were used to monitor a well during fracturing. The 
recordings of the geophones are then converted into arrival 
times and source location using an iterative, least square 

15 method. 

The present invention seeks to improve the amount of 
information gained from microseismic imaging of fracturing 
operations . 

20 

SUMMARY OF THE INVENTION 
The invention describes a waveform inversion algorithm to 
determine the mechanism, origin time and location of a 

25 seismic source. The algorithm is fully automatic, i.e. it 

does not require detection of seismic phases (such as P or S 
waves) or other parameters derived from data (e.g. 
polarization angles) . The algorithm is suitable for 
inversion in an arbitrary heterogeneous medium and takes 

30 advantage of a good velocity and density model, if it is 

available. An alternative version of the inversion algorithm 
(with location or origin time of the seismic source 
determined independently) can be used to invert for the 
mechanism of the source only. The algorithm uses reciprocity 

35 of the source and receivers by evaluating Green* s functions 
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in an arbitrary heterogeneous medium from the receiver 
locations. These Green's functions are then inverted to 
evaluate synthetic seismograms due to an arbitrary source 
mechanism from source locations. 

5 

Using preferably a grid search over all possible source 
locations and origin times, the full waveform synthetic 
seismograms are fitted to the data by the least-square 
method. The initial estimate of the origin time is set 
10 through cross-correlation of data and synthetics due to an 
arbitrary source mechanism. 

The inverted origin time is determined by a grid search 
around this initial estimate. The algorithm is robust to 
15 white noise added to the synthetic seismograms and is robust 
and particularly suitable for low frequency data in the 
frequency band from 0 Hz to 100 Hz, more preferably 0 Hz to 
50 Hz. 

2 0 These and further aspects of the invention are described in 
detail in the following examples and accompanying drawings. 

BRIEF DESCRIPTION OF THE DRAWINGS 

25 The invention will now be described, by way of example only, 
with reference to the accompanying drawings, of which: 

FIG. 1 shows a schematic illustration of fracturing 
operation; 

30 

FIG. 2 is a flowchart of steps performed in an example of 

the present invention; and ......... . r , m 
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FIG. 3 is a comparison of synthetic data with data derived 
using an example of the present invention. 
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(54) Abstract Title: Microseismic determination of location and origin time of a fracture generated by fracturing 
operation in a hydrocarbon well 

(57) In a method of monitoring a hydraulic fracturing 
operation in a hydrocarbon well 11, microseismic 
signals 14 originating from an induced fracture 14 
are recorded by three component geophone arrays 
121,131 in neighbouring wells 12,13. A wavefield 
inversion algorithm is used to determine the 
mechanism, origin time and location of the seismic 
source 14. The algorithm does not require the 
signals to be resolved into P-wave and S-wave data 
before inversion. The algorithm evaluates Green's 
functions, uses an existing velocity model and 
decomposes a moment tensor to yield parameters 
characteristic of the fracture. The recorded signals 
may be bandpass limited to 0-50 Hz or 0-100 Hz. 
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Method for Monitoring Hydraulic Stimulation 



This invention relates to methods for acquiring seismic data 
monitoring hydraulic stimulation such as fracturing rock 
5 layers to improve hydrocarbon production of a well. More 
specifically it relates to such methods using seismic 
methods to detect the location of an induced fracture. 



BACKGROUND OF THE INVENTION 
10 In certain situations, workers in the oil and gas industry 
perform a procedure known as "hydraulic fracturing" during 
oil exploration and drilling operations. For example, in 
formations where oil or gas cannot be easily or economically 
extracted from the earth, a hydraulic fracturing operation 
15 is commonly performed. Such a hydraulic fracturing operation 
includes pumping in large amounts of fluid to induce cracks 
in the earth, thereby creating pathways via which the oil 
and gas may flow. After a crack is generated, sand or some 
other material is commonly added to the crack, so that when 
20 the earth closes back up after the pressure is released, the 
sand helps to keep the earth parted. The sand then provides 
a conductive pathway for the oil and gas to flow from the 
newly formed fracture 

25 However, the hydraulic fracturing process does not always 

work very well. The reasons for this are relatively unknown. 
In addition, the hydraulic fractures cannot be readily 
observed, since they are typically thousands of feet below 
the surface of the earth. Therefore, members of the oil and 

30 gas industry have sought diagnostic methods to tell where 
the fractures are, how big the fractures are, how far they 
go and how high they grow. Thus, a diagnostic apparatus and 
method for measuring the hydraulic fracture and the rock 
deformation around the fracture are needed. 
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In previous attempts to solve this problem, certain methods 
have been developed for mapping fractures. For example, one 
of these methods involves seismic sensing. In such a seismic 
sensing operation, micro-earthquakes generated by the 
5 fracturing are analyzed by seismic meters, for example, 
accelerometers . 



A recent study on the use of microseismic imaging for 
fracture stimulation was published by J. T. Rutledge and W. 

10 S. Phillips. In an typical operational setting as described 
in greater detail in FIG. 1 below, three-component geophones 
were used to monitor a well during fracturing. The 
recordings of the geophones are then converted into arrival 
times and source location using an iterative, least square 

15 method. 

The present invention seeks to improve the amount of 
information gained from microseismic imaging of fracturing 
operations . 



SUMMARY OF THE INVENTION 
The invention describes a waveform inversion algorithm to 
determine the mechanism, origin time and location of a 

25 seismic source. The algorithm is fully automatic, i.e. it 

does not require detection of seismic phases (such as P or S 
waves) or other parameters derived from data (e.g. 
polarization angles) . The algorithm is suitable for 
inversion in an arbitrary heterogeneous medium and takes 

30 advantage of a good velocity and density model, if it is 

available. An alternative version of the inversion algorithm 
(with location or origin time of the seismic source 
determined independently) can be used to invert for the 
mechanism of the source only. The algorithm uses reciprocity 

35 of the source and receivers by evaluating Green's functions 
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in an arbitrary heterogeneous medium from the receiver 
locations. These Green's functions are then inverted to 
evaluate synthetic seismograms due to an arbitrary source 
mechanism from source locations. 



Using preferably a grid search over all possible source 
locations and origin times, the full waveform synthetic 
seismograms are fitted to the data by the least-square 
method. The initial estimate of the origin time is set 
10 through cross-correlation of data and synthetics due to an 
arbitrary source mechanism. 



The inverted origin time is determined by a grid search 
around this initial estimate. The algorithm is robust to 
15 white noise added to the synthetic seismograms and is robust 
and particularly suitable for low frequency data in the 
frequency band from 0 Hz to 100 Hz, more preferably 0 Hz to 
50 Hz. 



2 0 These and further aspects of the invention are described in 
detail in the following examples and accompanying drawings. 



BRIEF DESCRIPTION OF THE DRAWINGS 



25 The invention will now be described, by way of example only, 
with reference to the accompanying drawings, of which: 

FIG. 1 shows a schematic illustration of fracturing 
operation; 

30 

FIG. 2 is a flowchart of steps performed in an example of 
the present invention; and 



FIG. 3 is a comparison of synthetic data with data derived 
35 using an example of the present invention. 
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DETAILED DESCRIPTION 
A typical operational setting for monitoring hydraulic 
5 fracturing is illustrated in FIG 1 with a treatment well 11 
and geophone arrays 121, 131 located in neighboring wells or 
holes 12, 13. During the fracturing operation a fluid is 
pumped from the surface 10 into the well 11 causing the 
surrounding formation in a hydrocarbon bearing layer 101 to 
10 fracture. Acoustic waves 14 generated by the fracture 111 
propagate through the earth and are recorded by the three- 
components geophones of the two arrays 121, 131. 

For the present invention it is assumed that three 
15 components of the time history of particle velocity (or 

particle displacement) at several (N__r) downhole receivers 
were recorded during an acoustic emission. Furthermore, it 
is assumed the existence of an velocity model (of arbitrary 
complexity) of the volume of earth through which the seismic 
20 waves travels. The quality of the velocity model can be 

characterized by the length of time interval T_i (i= l..N_r) 
for which one is confident a synthetic seismograms can fit 
the data. These time intervals preferably include at least 
the S-wave arrival at all of the receivers. The use the 
25 particle displacement is preferred as it stabilizes the 

inversion as the particle velocity is more oscillatory than 
particle displacement. 

To find the relevant source parameters such as location 
30 vector x_s, origin time t_0 and moment tensor M, the misfit 
between a synthetic seismograms and data is minimized. In 
this inversion the misfit is defined by equation [1] : 



tl] 
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A = ££/ (d i (4,t-to)-V j (x,,4,t,M)) dt 

where d_j denotes a component of the particle velocity- 
recorded at the i-th receiver and U_j is the j-th component 
of the synthetic seismogram at the i-th receiver due to a 
5 source located at x__s characterized by a moment tensor M. To 
facilitate the description characters following a underscore 
appear as subscript in the equations . 

The source parameters that minimize equation [1] comprise 
10 the inverted solution. The j-th component of a synthetic 

seismogram at i-th receiver x x _r due to sources at locations 
x_s can be evaluated from the well known relation 

[2] 

L/}(x,,x*,t, M) = ^G fcjjm (x 5j x^)*M fcm (x,,^ 

15 

Here **" is a convolution in time, G_kj,m is the derivative 
of the Green's function along m-th coordinate axis and M_jk 
is a moment tensor of a point source located at x_s. 

20 The least-square minimum of the misfit given by equation [1] 
is in general non-unique. To alleviate this problem, it is 
preferred to make two assumptions: Firstly, approximating 
the source as a single point source x_s so that the sum over 
x_s in equation [2] disappears. Secondly, the source- time 

25 function can be approximated as a delta source- time function 
so that the convolution in the equation [2] is replaced by a 
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multiplication. Using these approximations the equation [2] 
reduces to 



[3] 

Uj(x s ,it*.,t } M) = G^x^x*,*) • M jk (x 3 ) 

+ G 2j}2 (x 3i x* , t) • M 22 (x s ) 

+ C?3j,3(x 5 ,x*,t) -M 33 (x 5 ) 

+ (C? 2 j,i(x 5 ,x*,t) +G ? i i , 2 (x 5 ,x*,t)) -M 21 (x 5 ) 

+ (<2 3 j,i(x,,x*,*) + G 3 j,i(x 5 ,x;^)) • M 3i (x 5 ) 

+ (^^(x^x*,*) + G Zji2 (xs,x). 7 t)) • M 32 (x 5 ). 

5 

It is known that equation [3] has a unique solution for M 
with a fixed origin time t__0, point-source location x_s and 
inversion model. Therefore, the trade-off among the source 
parameters can be minimized by a grid search over source 

10 locations and origin times for the best fitting moment 

tensors. The grid search for all possible origin times is 
numerically expensive and is therefore accelerated by 
estimating the origin time from cross-correlation of the 
synthetics and data and then using the grid- searching around 

15 this initial guess. The method used includes the following 
steps as illustrated in FIG. 2: 

- Following a recording of acoustic data from a fracture 
(Step 20) ; 

- estimate the initial origin time t_0(x_s) at 
20 every possible source location x_s (Step 21) ; 

- carry out a grid search around the estimated origin time 
for each source location (Step 22) . For each origin time 
find the unique solution M(x_s, t_0 (x__s) ) (least-square 
minimum) (Step 23) and evaluate the least-square misfit 

25 between the data and the synthetics (Step 24) ; and 

- store the best fitting solution for each source location 
(Step 25) . 
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The moment tensor of fracture together with the origin time 
and location can then be further evaluated (Step 26) as 
described below to find characteristics of the fracture. 

5 The initial estimate of the origin time is evaluated by 

cross-correlation of the data and synthetic seismograms for 
an chosen source mechanism, e.g. vertical strike-slip. The 
cross-correlation is evaluated over the time interval $(0, 
T_j) for each receiver j . The absolute values of the 

10 corresponding components for each receiver are cross- 
correlated and the time shifts of the maximum cross- 
correlation for each component are calculated. Using the 
absolute values of the seismograms for the cross-correlation 
reduces the dependency on the unknown source mechanism. The 

15 time shifts of each component and the known origin times of 
synthetic seismograms enables an estimation of the absolute 
origin time t°_i j for each component i and receiver j . The 
estimates are weighted by the maximum amplitude of the 
recorded seismograms to reduce poor estimates resulting from 

20 cross-correlating traces dominated by noise. It is worth 
noting that using the maximum amplitude as a weight in 
averaging the origin time assumes that the signal-to-noise 
ratio is proportional to the maximum amplitude of the 
recorded seismograms. The final estimate of the origin time 

25 is therefore an arithmetic weighted-average with weights of 
maximum amplitudes A_ij of i-th component at j-th 
receiver: 



t4] 




30 



57.0612 GB NP : I 

• ' • * • • ■ 



8 

This cross-correlation can be further improved at the 
expense of a more time intensive calculation by using the 
signal envelopes instead of the amplitudes. 



5 The true origin time is then found by grid- search around the 
initial estimate of the origin time within the dominant 
[shortest] period in the signal. The limiting of the grid 
search to the dominant period of the signal requires the 
initial estimate of the origin time [4] to be within the 

10 dominant period. This is typically the case for the S-wave 
arrival. The grid search around the initial estimate of the 
origin time [4] eliminates the problems with the cycle- 
skipping as the cross-correlation function tends to peak 
every 1/2 -period of the dominant period (usually the minimum 

15 period present in the data) . 



The length of the time step in the grid search is set to 
obtain the required accuracy of the misfit [1] . Assuming 
that the synthetic seismograms match the data (i.e. using 
20 the true moment mechanism and evaluating the synthetic 
seismograms in the true model from the true source 
location) , normalized misfit of a harmonic signal with 
period T, due to a time shift of OCT in the origin time, can 
be evaluated as 



25 



[5] 



E = /o r NnM)-sin(.(, + aT))]^ = ± _ 
2f*[sm(uji)] 2 dt 



30 



The definition of error in equation [5] has a maximum of 2 
for 1/2 period shift and even a small time shift causes a 
large error for a misfit defined analogously to equation 
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[1] . The length of time step for the grid search can be set 
to 2otT for which the maximum error of evaluation of misfit 
reaches a certain limit. For example, a shift of 0.05 T$ 
(a=0.05) may cause relative error E=0.05. Thus, a search for 
5 origin time with a grid step of 0.1T (T is the dominant 
period in my seismograms) should not cause an error of 
evaluation in the misfit function larger than 0.05. 



The last part of the method is to identify a unique solution 
10 M(x_s, t_0(x_s)) for each origin time and source location. 
It is known that the moment tensor with the least-square 
minimum fit of the equation [1] is: 

[6] 



Mi(x a ) - (A-%(x,)Dj(ic s ). 

15 

Here M_I(bar) is the i-th component of six elements vector: 
M_I(bar) = M_ll, M_2 (bar) = M_12= M_21, M_3 (bar) = M_22, 
M_4(bar) = M_13 = M_31, M_5(bar) = M_3 = M_32, M_6 (bar) = 
M_33, and D has six independent elements 

20 

[7] 

N r 3 „Ti 

i=0 j=0 J ° 

Here k=0 . . . 5 and g_jk is defined by the following notation: 



25 



[8] 



57.0612 GB NP : .' .\ . 

■ » t i tit • 

• it t t 

• • • f t 

10 

9jl{ x sj x r j t) — Gy^ifXj, Xr, t) 

^2(^5, X rj t) = G 2 j,l( x s, x r>£) + Gij 3 2(x s ,x r7 ^) 

9js{ x si Xr, t) = G f 2j,2(x 5 , X r , t) 

g j4 (x s ,Xr,i) = G 3 j,i(x^x r ^) +Gy j3 (x„x rj t) 

fly"6 ( x s > x r > = ^37,3 ( x s > x r ? • 



Finally, A is a 6x6 matrix with elements: 



[9] 

N r 3 . Ti 

4w(x,) = / 

i=Q j=0 J ® 

5 

The integration steps of [7] and [9] can be accelerated by 
using a time window t_min to t_max, where t_min is a time of 
arrival of a first energy from the source ( fracture) as 
10 identified by an event detector and t_max is the maximum 

time for which the waveforms are matched, e.g., the time of 
arrival of the phase with maximum amplitude. This 
modification excludes the effect of reflections or tube 
waves in the recorded data. 

15 

It is further feasible to associated with every recording 
device or trace a weighting function that indicates the 
quality of the receiver and/or recorded data. These weights 
could be introduced into the present equations [7] and [9] . 



The synthetic Green's function in equation [3] is then 
evaluated by computing three times N_r full waveform 
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simulations (using a finite-differences) . For each three- 
component receiver, three responses due to three orthogonal 
single force sources at the receiver positions are computed 
and derivatives of the velocity (or displacement) are stored 
5 at every possible source location, x_s . The synthetic 

seismograms are evaluated with a delta function as a source- 
time function. Using reciprocity, derivatives of Green's 
functions for every possible source location to every 
receiver position are evaluated. Equation [3] shows that six 
10 traces at every possible source location must be stored. 

The above equation provides a complete set of steps to 
calculated the moment tensor M from three component 
recordings of the wavefield. The tensor itself is then 

15 decomposed to yield parameters characteristic of the 

fracture. Methods to decompose the moment tensor M have been 
developed for the purpose of analyzing earthquakes and are 
described for example by V. Vavrycuk in: Journal of 
Geophysical Research, Vol 106, No B8, August 10, 2001, 

20 16,339-16,355. The parameters obtained by such decomposition 
include the normal of the fracture n, the slip direction N, 
and products of the Lame coefficients with the slip u of the 
fracture, i.e., Jiu and Xu respectively. Alternatively, the 
moment tensor can be inverted for a set of parameters 

25 including the orientation of the pressure P and tension T 

axes, parameter K = X/\i and inclination a of the slip u from 
the fracture. These parameters provide information on the 
fracture orientation and slip direction which in turn can be 
used to control the hydraulic fracturing operation. 



30 



The accuracy of the inversion from recorded data d_j to the 
moment tensor M of the source can be further improved by 
bandlimiting the frequency of the data. While restricting 
data to the 0-100 Hz band yields satisfactory results, an 
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improved accuracy is gained by limiting the data further to 
the 0-75 Hz and even the 0- 50 Hz band. In FIG. 3 there is 
shown a plot of (synthetic) geophone velocity measurements 
31 in x, y and z directions overlaid with the corresponding 
5 . traces 32 re-calculated using the moment tensor derived by 
the method described above (with a known velocity model) . 
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CLAIMS 

A microseismic method of monitoring fracturing 
operation in hydrocarbon wells comprising the steps of 
obtaining multi-component signal recordings from 
locations in the vicinity of a facture; and performing 
a wavefield inversion to determine parameters 
representing a source mechanism of said fracture. 

The method of claim 1 wherein the signal recordings are 
at least for the purpose of determining the moment 
tensor bandlimited to 0 to 100 Hz . 

The method of claim 1 wherein the signal recordings are 
at least for the purpose, of determining the source 
mechanism bandlimited to 0 to 5 0 Hz. 

The method of claim 1 wherein the signal recordings are 
at least for the purpose of determining the source 
mechanism bandlimited to 5 to 50 Hz. 

The method of claim 1 wherein the signal recordings are 
not resolved into P-wave or S-wave event prior to the 
wavefield inversion. 

The method of claim 1 wherein the parameters of the 
source mechanism are represented by a moment tensor. 

The method of claim 1 including the step of minimizing 
the difference between recorded signals and synthetic 
signals . 

The method of claim 1 including the step of minimizing 
the difference between recorded signals and synthetic 
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signals, wherein said synthetic signals depend on 
source location, receiver location, time and source 
mechanism, 

5 9. The method of claim 8 wherein the step of minimizing 
the difference between recorded signals and synthetic 
signals includes a grid search over source locations, 
origin times for an optimized source mechanism. 

10 10. The method of claim 1 further comprising the step of 
determining location and origin time of a source 
represent the fracture. 

11. The method of claim 10 further comprising the steps of 
15 - estimating the initial origin time at possible source 

locations ; 

- carrying out a grid search around the estimated 
origin time for each source location;. 

- for said origin times finding the unique solution of 
20 a moment tensor of the source; 

- evaluating the least-square misfit between the 
recorded signals and synthetic signals derived by 
calculating the signals caused by a source of said 
moment tensor at the receiver locations; and 

25 - storing the best fitting solution for each source 

location. 
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12. The method of claim 1 wherein a source time function of 
the fracture is approximated by a delta function. 

13. The method of claim 1 including the step of evaluating 
a Green's function to derive the source mechanism from 
the recordings. 
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